load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


maskFileList = systemfunc("ls /public/users/wrf/Data/SanXia/*_mask_d02.csv")
nRegion = dimsizes(maskFileList)

if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/sanxia_110820/wrfout_d02_2011-08-20_00:00:00"
    print("input_file is undefined, using default: " + input_file)
end if

if (.not. isvar("output_prefix")) then
    print("output_prefix is undefined")
    output_prefix = "~/wrfpost/test/sanxia_precip"
    print("output_prefix is undefined, using default: " + output_prefix)
    system("mkdir " + output_prefix)
end if

; default precipitation delta-t is 1 hour
if (.not. isvar("precip_dt")) then
    precip_dt = 6 ; output 1h precip
    print("precip_dt is undefined, using default: " + precip_dt)
end if


output_header = "Time"

; Read areaMaskFile
do i_region = 0, nRegion-1, 1
    areaMaskFile = maskFileList(i_region)
    region_name = readAsciiHead(areaMaskFile, 1)
    wrf_mask = readAsciiTable(areaMaskFile,4,"float",1)
    ij_masks = tointeger(wrf_mask(:,2:3))
    dim_mask = dimsizes(ij_masks)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;           Open NetCDF File               ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    a = addfile(input_file + ".nc", "r")
    filename = systemfunc("basename " + input_file)

    domain = substring(filename, 7, 9)
    start_date = substring(filename, 11, 23)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;              Process Time                ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    times = a->Times
    hours = WRF_Times2Udunits_c(times, 0)
    ntimes = dimsizes(hours)
    wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)

    ; UTC -> CST
    hours = hours + 8 ; 
    timeStr_cst = ut_string(hours, "%Y-%N-%D %H:%M:%S")

    dt = round(precip_dt / wrf_dt, 3)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;            Read Precipitation            ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    rain_exp = wrf_user_getvar(a,"RAINNC",-1)
    rain_con = wrf_user_getvar(a,"RAINC",-1)
    rain_tot = rain_exp + rain_con
    dim_rain = dimsizes(rain_tot)

    rain_crop = new((/dim_rain(0),dim_mask(0)/), "float")
    do i = 0,dim_mask(0)-1,1
        rain_crop(:,i) = rain_tot(:,ij_masks(i,1), ij_masks(i,0))
    end do
    rain_crop_avg = dim_sum(rain_crop) / tofloat(dim_mask(0))
;    rain_tend = array_append_record(0.0, (rain_crop_avg(dt:ntimes-1:dt) - rain_crop_avg(0:ntimes-1-dt:dt)), 0)
;    timeStr_tend = timeStr_cst(0:ntimes-1:dt)
    rain_tend = rain_crop_avg(dt:ntimes-1:dt) - rain_crop_avg(0:ntimes-1-dt:dt)
    timeStr_tend = timeStr_cst(dt:ntimes-1:dt)
    output_header = output_header + ", " + region_name
    if (.not. isvar("output_body")) then
        output_body = timeStr_tend + ", " + sprintf("%6.3f", rain_tend)
    else
        output_body = output_body + ", " + sprintf("%6.3f", rain_tend)
    end if
    delete(wrf_mask)
    delete(ij_masks)
    delete(rain_crop)
end do

output_text = array_append_record(output_header, output_body, 0)
output_filename = output_prefix + "/sanxia_areal_precip_" + domain + "_" + start_date + ".csv"
asciiwrite(output_filename, output_text)
